#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

using namespace std;

// Random Number generation with Maxwell distribution
double FOverG(double y){
 return y*exp(-0.5*y*y + y);
}

double randMaxwellDist(gsl_rng *r){
  double Y = 0.0;
  double U = 0.0;
  double c = 2.2039;
  double X = 0.0;

  do {
    Y = gsl_ran_exponential(r,1);
    U = gsl_rng_uniform(r);
  } while(U > FOverG(Y)/c);

  X = Y;
  return X;
}


void randMaxwellDist2(){
 const gsl_rng_type *T; T= gsl_rng_default;
 gsl_rng *r; r = gsl_rng_alloc(T);
 long unsigned int seed;
 seed = 100;
 int N = 1000;
 gsl_rng_set(r,seed);
// for (int i=0; i<N;i++)
 int count = 0;
 double Y = 0.0;
 double U = 0.0;
 double c = 2.2039;
 double X = 0.0;
 while (count < N){
  Y = gsl_ran_exponential(r,1);
  U = gsl_rng_uniform(r);
  if ( U <= FOverG(Y)/c ){
    X = Y;
    cout<< X<<endl;
  }
 // cout<<0.01*count<<" "<<FOverG(0.01*count)<<endl;
  count++;
 }
 gsl_rng_free(r);
}
/*
double *function(int const N){
  double *array;
  array = new double[N];
  for (int i=0;i<N;i++){
    array[i] = 0.1*i-0.5
  }
  return array;
}
*/
int main(){
  const gsl_rng_type *T; T = gsl_rng_default;
  gsl_rng *r; r = gsl_rng_alloc(T);
  long unsigned int seed;
  seed = 100;//(long unsigned int) time(NULL);
  gsl_rng_set(r,seed);
  int N = 10000000;
  for (int i=0; i<N; i++)
    cout<<randMaxwellDist(r)<<endl;
  
}


